To appear in the Springe r-Verlag IMA volume on Compatible Discretizations , 
Eds. Arnold, Bochev, and Shashkov, 2005 

ON THE ROLE OF INVOLUTIONS IN THE DISCONTINUOUS 
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MAGNETOHYDRODYNAMIC SYSTEMS 

TIMOTHY BARTH* 

Abstract. The role of involutions in energy stability of the discontinuous Gaierkin (DG) dis- 
cretization of Maxwell and magnetohydrodynamic (MHD) systems is examined. Important differ- 
ences are identified in the symmetrization of the Maxwell and MHD systems that impact the construc- 
tion of energy stable discretizations using the DG method. Specifically* general sufficient conditions 
to be imposed on the DG numerical flux and approximation space are given so that energy stability is 
retained. These sufficient conditions reveal the favorable energy consequence of imposing continuity 
in the normal component of the magnetic induction field at interelement boundaries for MHD dis- 
cretizations. Counterintuitively, this condition is not required for stability of Maxwell discretizations 
using the discontinuous Gaierkin method. 
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1. Overview. Various mathematical models such the Maxwell equations 
governing electrodynamics and the magnetohydrodynamic (MHD) equations 
modeling fluid plasmas have the added complexity of possessing involutions . An 
involution in the sense of conservation law systems is an additional equation that 
if satisfied at some initial time is satisfied for all future time for both classical 
and weak solutions [Boi88, Daf86]. Involutions should not be confused with con- 
straints that are needed for closure of the system. An example of such a constraint 
is the continuity equation in incompressible flow. In this note, the role of involu- 
tions in obtaining energy stable discretizations using the discontinuous Gaierkin 
method [RH73, LR74, JP86, CLS89, CHS90] is briefly examined. Specifically, 
the surprisingly different role played by involutions in the discontinuous Gaierkin 
(DG) discretization of Maxwell and ideal compressible MHD systems is con- 
trasted. Although both systems possess solenoidal involutions, it is the interplay 
between involutions and symmetrization of the Maxwell and MHD systems that 
enters fundamentally into the construction of stable discretizations. In this regard, 
the two systems are vastly different. The Maxwell equations are naturally ex- 
pressed in essentially symmetric form. Consequently, the analysis given in Sects. 
2. 1 and 3. 1 shows that “standard” DG discretizations can then be used. In contrast, 
symmetrization of the MHD system utilizes the solenoidal involution as a nec- 
essary ingredient in the symmetrization process. Details of this symmetrization 
process are given in Sect. 2.2. Thus, the precise sense in which involutions are sat- 
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isfied in element interiors and across interelement boundaries enters prominently 
into the MHD discrete energy analysis. The analysis of Sect. 3.2 gives general 
sufficient conditions to be imposed on the DG numerical flux and approximation 
space in the presence of involutions so that energy stability is retained. These 
sufficient conditions reveal the favorable consequences of imposing continuity in 
the normal component of the magnetic induction field at interelement boundaries 
for MHD discretizations. This is a condition that is not required for stability of 
Maxwell discretizations using the discontinuous Galerkin method but is often a 
requirement of other methods that build satisfaction of solenoidal conditions into 
the discretization. Techniques for achieving this include staggered mesh and spe- 
cialized differencing techniques [Yee66] as well as edge, face, and volume finite 
element formulations [Ned80, Bos98, BR02] or the discrete mimetic approxima- 
tions as given in [HS99]. The present analysis for MHD also provides alternatives 
to the “divergence cleaning” procedures designed to exactly or approximately sat- 
isfy the solenoidal condition, see [BB80, TOO, DKK + 02, BK04] and references 
therein. Since the DG method reduces to the simplest finite volume method in 
the special case of piecewise constant basis approximation, the results given here 
impact finite volume discretization as well. 

2. Symmetrization of Conservation Laws without Involution. Consider 
the Cauchy initial value problem for a system of m coupled first-order differential 
equations in d space coordinates and time which represents a conservation law 
process. Let u(x 5 t) :JR d x IR + »-> IR m denote the dependent solution variables 
and f(u) : IR m »-»• lR Tnxd the flux vector. The model Cauchy problem is then 
given by 

n u f + ft = 0 

K } \ u(x, 0) = uoOr) 

with implied summation on the index i ~ 1, . . . , d. Additionally, the system is 
assumed to possess a convex scalar entropy extension. Let U( u) : IR m »-+ IR and 
F( u) : lR m i ♦ ]R d denote an entropy-entropy flux pair for the system such that 
in addition to (2.1) the following inequality holds 

(2.2) U t + Fi, Xi < 0 

with equality for classical (smooth) solutions. In the symmetrization theory for 
first-order conservation laws without involution [God61, Moc80], one seeks a 
mapping u(v) : ]R m i-* IR™ applied to (2.1) so that when transformed 

(2.3) u )V v )t + f ifV v = 0 

the matrix u >v is symmetric positive definite (SPD) and the matrices f ijV are 
symmetric. Clearly, if twice differentiable functions U(v) : 3R m i-> IR and 
Fi(v) : IR m i— ► IR can be found so that 

(2-4) u = W£, f< = Tj y 
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then the matrices 


U v — ^,vvj fi,v — i 


t,vv 


are symmetric. Further, we shall require that U{\) be a convex function such that 
(2.5) 


,i» m = + ~ 


v * oo | v I 

so that U (u) can be interpreted as a Legendre transform of U(v) 
U(u) = sup{v ■ u - U{v)} . 


From (2.5), it follows that 3 v* € IR m such that v-u-W(v) achieves a maximum 
at v* 

(2.6) (7(u) = v* • u — ZY(v*) . 

At this maximum u = V (v*) which can be locally inverted to' the form v* — 

v(u). Elimination of v* in (2.6) yields the simplified duality relationship 

U (u) = v(u) • u - W(v(u)) . 


Differentiation of this expression 

(2.7) £/„ = v 4- v jU u - v, u ^ = v 

gives an explicit formula for the entropy variables v in terms of derivatives of the 
entropy function U( u). Using the mapping relation v(u), a duality pairing for 
entropy flux components is defined 

Fi( u) = v(u) • fi(u) - ^i(v(u)) . 

Differentiation then yields the flux relation 

F it u = V • f iiU + V, u fi - V,nF^ v = V • f <iU 

and the fundamental relationship for classical solutions 

v • (u t + = u\ t + Fi jXi — 0 . 

These relationships are used extensively in the discrete energy analysis of the 
discontinuous Galerkin method. 

2.1. Maxwell Equations in Symmetric Form. The time-dependent 
Maxwell equations are given by 

H (b) + v x ( _ E B ) = ( Maxwe11 equations) 
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where E e IR d , B G IR d , p c G H, and j G IR d denote the electric field, magnetic 
induction, charge and current density with eo anc * c the free-space permittivity and 
speed of light, respectively. If the charge conservation equation 

(2.8) (pc),t + V * j = 0 

is satisfied for all time then the Maxwell system possesses the following involu- 
tions 


V • E = p c /e o 

V • B = 0 . 


Writing the Maxwell system in matrix coefficient form reveals that the above 
system is essentially already in symmetric form using the variables u = (E, B) r 


u )t + Ai u x . =q(u) , 



c 2 Mi 

0 


where in three space dimensions 


"0 

0 

o' 


*0 

0 

-1" 


" 0 

1 

0“ 

0 

0 

1 

, M 2 = 

0 

0 

0 

, M 3 = 

-1 

0 

0 

0 

-1 

0 


1 

0 

0 


0 

0 

0 


Consequently, a suitable entropy -entropy flux pair for the Maxwell system are 
given by the scaled “square entropy” and square entropy flux 


u( u) = i(|E | 2 + c 2 |B| 2 ) , F( u) = c 2 (E X B) . 


Using this entropy function, the symmetrization variables and right symmetrizer 
are then obtained 


( E \ 


Idxd 

U 2 b) 

i u >v — 

C~ 2 Idxd_ 


thus rendering the coefficient matrices symmetric as expected 


AiU v = 


0 

M[ 


Mi 

0 


Observe that the Maxwell system has been successfully symmetrized without uti- 
lizing the involutions. Consequently, the energy analysis for Maxwell’s equations 
in a vacuum domain is identical to the energy analysis for conservation law sys- 
tems without involution as also observed in [CLS04]. 
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2.2. Ideal MHD in Symmetric Form. The equations of ideal compressible 
MHD are given by 


d_ 

dt 


( p \ 

pV 




+ V- 


where p 6 IR, V 6 
magnetic induction. 


( pV \ 

pVV + I dxd (p + |B| 2 /2) — B B 
(J5 + p + |B| 2 /2) V — (V • B) B 

VB - BY J 


(Ideal MHD) 


lR <i , B € lR d , and p € IR denote the fluid density, velocity. 


and pressure with E G IR the total specific energy given by 


E = ^- 1 + p \ V \ 2 /2 + \ B \ 2 /2 


and 7 the ratio of specific heats. In addition, the MHD system possesses the 
solenoidal involution 


V * B = 0 

which in consistent with the absence of experimentally observed magnetic 
monopoles. 

It is well known that thermodynamic entropy s is transported along velocity 
induced particle paths for ideal MHD. Recall that 5 = log (pp -7 ) for MHD so 
that a differential of 5 is given by 

7 1 

ds = — dp H — dp . 

p p 

Inserting equations derived from the MHD system (2.2) yields 

3 t + V ■ V s + (7 - 1 ) ^ V • B = 0 

P 

or after combining with the continuity equation 

(ps),t + div(pVs) + (7 - 1 )-- V B V • B = 0 

P 

suggesting that U (u) = - ps may be a suitable entropy function only if the in- 
volution V • B = 0 is satisfied. Indeed, a straightforward calculation for ideal 
MHD shows that this entropy function does not symmetrize the system under the 
change of variable u i—> v with v = (see for example Barth [Bar98]) 

f V * f T y 

since the involution equation has not been used. Godunov [God72] observed this 
phenomenon as well which lead to his development of a symmetrization tech- 
nique for ideal MHD. The basic technique is reviewed here using a modified pre- 
sentation from that originally given. The model MHD system with solenoidal 
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involution is given by 


(2.9) 


u ,t + ft, it “ 0 

B i jXi = 0 
u(x,0) = u 0 (x) 


with convex entropy extension 

(2.10) U ti + Fi, Xi < 0 . 

To analyze this system, Godunov considered augmenting the MHD system by 
adding multiples of the involution where the multipliers are themselves the gra- 
dient of a scalar function <£(v) : IR m IR with respect to the symmetrization 
variables v 

u ,t + fi,ii + Bi )Xi = 0 . 

Consider the following ansatz for the dependent variables u and flux components 
ft 


u = U T v 

- r ( v ) Bi 

with U a convex scalar function and r(v) : IR™ i-> IR m an unknown vector- 
valued function. Observe that the augmented MHD system 

(2.11) (Wv), t + (^. v -r(v)B 1 )f Ji +^ v B i , li =0 

possesses a symmetric quasilinear form in v variables whenever r(v) = (f) T v since 
the system (2.1 1) then reduces to 

^,vv "T («^i,vv “ ^,vvBi) V (Ii = 0 

SPD SYMM 

so that the final flux relationship is obtained 

The entropy function U(u) for MHD can be interpreted as a Legendre transform 
of U(v) 


JJ{ u) = sup {v • u — £Y(v)} 

v 


eventually producing the generalized duality relationships 

U (u) = v(u) • W,v(v(u)) - W(v(u)) 

(2.12) Fi{u) = v(u) • Jf ijV (v(u)) - Ji(v(u)) 
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so that for classical MHD solutions 

v ■ (u t + + (j) T v B i>Xi ) = U >t + F i}Xi = 0 . 

This relationship will be used heavily in later analysis of the discontinuous 
Galerkin method. 

Choosing the entropy function (7 (u) = —ps yields 0(v) = ( 7 — l)pV-B/p, 
a homogeneous function of degree one in v so that <p = v. The resulting in- 
volution multipliers 4> y are identical to those derived by Powell [Pow94] using 
a completely different argument motivated by (in part) the lack of Galilean in- 
variance of the original MHD system and the subsequent addition of a divergence 
wave family into the local Riemann problem solution to restore Galilean invari- 
ance. 

REMARK 2.1. Observe that MHD provides one particular example of a 
symmetrizable system with a given entropy- entropy flux pair {U, F x } for which 
the flux is not expressed as the gradient of a primative function Ti but rather 

U = Tfy - 4>l . 

In fact, for the specific MHD entropy function U (u) — — ps, it is possible to show 
that there cannot exist a function T x such that 


Thus , the DG energy analysis of MHD systems is fundamentally different from the 
energy analysis of systems not possessing involutions. 

3. The DG Finite Element Method. Let ft denote a spatial domain com- 
posed of stationary nonoverlapping elements K x , ft = U K x , K X C\K^ = 0 , i ^ j 
and time slab intervals I n = t" +1 ], n = 0, . . . , N — 1. Both continuous in 

time approximation and full space-time approximation on tensor space-time ele- 
ments Ki x I n will be considered in the analysis. It is useful to also define the 
element set T — {Ki , . . .} and the interface set 8 = {ei , e 2 5 • ■ •} with inter- 

face members K x r\ Kjfl ^ j of measure d — 1 corresponding to edges in 2-D 
and faces in 3-D. Let Vk{Q) denote the set of polynomials of degree at most k 
in a domain Q C IR d . In the discontinuous Galerkin method, the approximating 
functions are discontinuous polynomials in both space and time 

\> h = {w | W| Kx/ „ € ( VkiKxnf ,\/KeT,n = 0,...,N-l} . 

Alternatively, [CLS89, CHS90, Shu99] utilize a semi-discrete formulation of the 
DG method together with Runge-Kutta time integration. In this case, the set of 
approximating functions are discontinuous polynomials in space and continuous 
functions in time denoted by V 

For ease of exposition, the spatial domain ft is assumed either periodic in all 
space dimensions or nonperiodic with compactly supported initial data. In this 
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domain, we first consider the standard first-order Cauchy initial value problem 
(without involution) 


(3.1) 


( u ,t + U,xi — o 

\ U = U 0 (x) 


with convex entropy extension 


(3.2) U t + F lyXi < 0 . 

The DG method for the time interval with weakly imposed initial data 

t °_ ) obtained from a suitable projection of the initial data v(uq(x)) is given 
by the following statement: 

DG FEM: Find \ h € V h such that 

(3.3) B DG (vh,Wh.) = 0 , Vw k € V h 


with 


Bdg(v,w)=^I / / -(u(v)-w, t + fi(v)-w , Xi )dxdt 

n= 0 \K€T J K 

4- / / w(z_) * h(v(x_), v(z + ); n) ds dt 

£TirJl«JdK 


(3.4) 


KZT U1 *' an 

+ Y. f (w(t" +1 ) -u(v(C +1 )) - w(t") -u(v(t”))) dx j 
Ker^ K J 


with suitable modifications when source terms are present. In this statement 
h(v_,v + ;n) : IR m x IR m x lR d ]R m denotes a numerical flux function, 
a vector- valued function of two interface states v ± and an oriented interface nor- 
mal n with the following consistency and conservation properties: 

• Consistency with the true flux, h(v, v; n) = f (v) • n 

• Discrete cell conservation, h(v_, v + ; n) = — h(v +5 v_; — n) . 

For a given symmetrizable system with entropy function U( u), the DG method 
is uniquely specified once V h , the entropy function U (u), and the numerical flux 
function h(v_,v + ;n) are chosen. In this formulation, the finite-dimensional 
space of symmetrization variables are the basic unknowns in the trial space 
V h and the dependent variables are then derived via u(v/ l ). When not needed for 
clarity, this mapping is sometimes explicitly omitted, e.g. f/(v^) is written rather 
than U (u(v/ 1 )). An important product of the DG energy analysis given below are 
sufficient conditions to be imposed on the numerical flux so that discrete entropy 
inequalities and total entropy bounds of the following form are obtained for the 
discretization of the Cauchy initial value problem: 

• A local cell entropy inequality assuming continuous in time approxima- 
tion, v h e v£ 

[ U(vk)dx+ f F(v- i h^+y ) n) ds < 0 , for each KeT 
dt Jk JdK 

(3.5) 
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where v +)/l ; n) denotes a conservative numerical entropy flux. 

Summing over all elements then yields the global inequality 

(3.6) j t J^U(v h )dx< 0. 

• A total entropy bound assuming full space-time approximation, 6 V h 

f U(u*(t°_)) dx < f U(u(v k (x,t^))) dx < f C/(u(v /l (x,f° ))dx 
Jn Jo. Jo 

(3.7) 

where u*(f° ) denotes the minimum total entropy state of the projected 
initial data 


Under the assumption that the symmetrizer u v remains spectrally 
bounded in space-time, i.e. there exist positive constants c 0 and Co in- 
dependent of Vh, such that 

0 < c 0 ||z|| 2 < z • u i v(v /l (i,t))z < Co ||z|| 2 

for all z / 0, the following L 2 stability result is then readily obtained 
for the Cauchy problem 


||u(v ft (.,^))-u*(t 0 _)|U 2 ( n )< Hv h (-,t 0 _))-u*(t°_)|U 2(n) 

3.1. DG Energy Analysis for Systems without Involution. In this section, 
the DG energy analysis for systems of conservation laws without involution is 
reviewed. From Sect. 2.1 it was shown that this analysis is also the relevant 
analysis for the Maxwell system since this system can be symmetrized without 
using the Maxwell system involutions. Consequently, consider the DG method 
applied to the nonlinear system (3.1). For brevity, we avoid the introduction 
of trace operators and instead use the shorthand notation for interface quanti- 
ties f± = /(v(i±)), (f)t = (/_ + /+)/ 2 and [f}± =/+-/_. An energy 
analysis assuming continuous in time functions, G V f, yields the following 
cell-wise local entropy inequality which build upon previous scalar conservation 
law analysis for DG by [JJS95, JS94] and further related DG analysis for systems 
in [CS97] and [Bar98, Bar99]. 

Theorem 3.1 (DG Semi-Discrete Cell Entropy Inequality). Letv h e 
denote a numerical solution obtained using the discontinuous Galerkin method 
(3.4) assuming a continuous in time approximation for the Cauchy initial value 
problem (3.1) with convex entropy extension (3.2). Assume the numerical flux 
h(v_, v+; n) satisfies the system E-flux condition 

[v]t • (h(v_,v + ;n) -f(v(0)) • n) < 0 , V6> 6 [0, 1) 


(3.8) 
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where v(9) = v_ -f 0 [v]l. The numerical solution then satisfies the local 
semi-discrete cell entropy inequality 

(3.9) 4 [ U{v h )dx+ [ F(v_ jh ,v +ih ;n)ds < 0 , for each K E T 

dt Jk JdK 

with 

(3.10) F(v_,v+;n) = (v)± ■h(v_,v + ;n) - (F ■ n)t 
as well as the global semi-discrete entropy inequality 

(3.11) j t jjJ{w h )dx< 0. 


Proof \ Evaluate the energy, 5 dg(v^, v^), for a single stationary element K 
assuming continuous in time functions 


[ v ■ u £ dx = — [ U dx 

Jk ’ Jk 

-a 


“V iI{ • fi dx + 


v_ h ds 


— - ( / dx+ v_ • h < 

— - / (— • n 4- v_ ■ h) ds 

JdK 

= - [ ( F(v_,v + ;n) + D(v_,v + ;n) ) 

JdK ' v ' ' / 

Conservative Flux Entropy Dissipation 


ds 


for carefully chosen conservative entropy flux and entropy dissipation functions 

F(v_,v + ;n) = (v)t • h(v_, v + ; n) - {F-n)t. 

D(v_,v + ;n) = — i([v]± • h(v_,v + ;n) - • n]±) . 

Observe that the chosen form of F(v_,v+; n) is a consistent and conservative 
approximation to the true entropy flux F(v) 

• F(v, v; n) = (v • f - F) • n = F • n (consistency) 

• F(v_, v + ;n) = -F(v+,v_;-n) (conservation) . 

The only remaining task is to determine sufficient conditions in the design of the 
numerical flux h(v_,v + ;n) so that D(v_, v + ; n) > 0. Rewriting the jump term 
appearing in the entropy dissipation term as a path integration in state space 

£>(v_,v + ;n) = — l([v]l ■ h(v_, v+; n) — [F ■ n]t) 
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= -|[v]t • (h(v_,v+;n )~ J q ^( V W) • n ddj 

= -i[v]t • ^h(v_,v+;n) - f(v(0)) • n dB^j 

= ~\J Q M± •(h(v_,v + ;n)-f(v(0))-n) d6 . 

A sufficient condition for nonnegativity of D(v_,v + ;n) and the local cell en- 
tropy inequality (3.9) when applied to finite-dimensional subspaces is that the 
integrand be nonpositive. This yields a system generalization of Osher’s famous 
E-flux condition for scalar conservation laws given in [Osh84] 

(3.12) [v]± • (h(v_,v + ;n) -f(v(6>)) • n) < 0 , V0 € [0, 1] . 

Summation of (3.9) over all elements in the mesh together with the conservative 
telescoping property of F(v_ , v + ; n) yields the global entropy inequality (3. 1 1). 
□ 

Let Aj < A 2 < ■ • • < A™ denote ordered eigenvalues of f u . Some specific 
examples of system E-fiuxes (proofs omitted here) include 

• Symmetric variable variant of the local Lax-Friedrichs flux 

(3.13) h SLF (v_,v + ;n) = (f • n)+ - ^A max [u(v)]*+ 

with 


A max = sup max |A £ (v(^))| 
o<£<i i<*<™ 


where v(£) = v_ + £ [v]t. 

• Symmetric variable variant of the Harten-Lax-van Leer-Einfeldt flux 
[HLvL83, EMRS92] 

(3.14) h S HLLE(v_,v + ;n) = {f • n ) ± - ^^ HlLE (v.,v + ;n) 
with 

h d r \ A max + A m in 2A max A m in , 

shlle( v — » v + ’ ) = T — r — [f(v;n)]r~- — — [u(v)]t 

■^max ''min ^max / 'min 

and 


Amax — sup max(0, A m (v(£))) , A min = inf min(0, Ai(v(f))) 
0<£<1 0<^<1 


where v(£) = v_ -f £ [v] 
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Fully discrete entropy bounds are readily derived assuming DG finite element 
discretization in time. 

Theorem 3.2 (DG Fully-discrete Total Entropy Bounds). Let v h e 
V h denote the space -time numerical solution obtained using the discontinuous 
Galerkin method (3.4) for the Cauchy initial value problem (3.1) with convex en- 
tropy extension (3.2). Assume the numerical flux h(v_ , v+ ; n) satisfies the system 
E-flux condition 

[v]l • (h(v_,v + ;n) - f(v(0)) • n) < 0 , V0 e [0, 1] 

where v(0) = v_ + 8 [v]j\ The numerical solution Vh then satisfies the total 
entropy bound 

(3.15) [ U(u*(t°_ )) dx < / U(u(v h (x,t^)))dx < f U(u(v h (x, $°)) dx 
Jo Jo Jo 

where u* (t°_) denotes the minimum total entropy state of the initial projected data 

u ‘ (< °-> = SS=(nj/ n u(v ‘ (a '’ ,0 - ))<it - 


Proof. Analysis of the spatial terms follows the same path taken in Theorem 
3.1 (omitted here) with an additional integration performed in the time coordi- 
nate. Consider the energy of the remaining time evolution terms in (3.4) after 
integration-by-parts for a single time slab interval I n 

[ f vu yt dxdt+ [ v{t+) • [u]*t dx— [ f Ujdtdx+ f v(t+) - [u}\idx 
Ji n Jo ’ Jo ~ JoJi n Jo 

= [ 

Jo. 

Taylor series with integral remainder together with the duality relationship (2) 
yields 




[tn?-v(e) 


[u]ji+£ n = 0 , R n = [\l- 

J o 


• e ) 


r(v(*)) 


de > o 


t n 

where v(0) = v(fl) + 0 [v] t « . Inserting into the time evolution terms 


[ [ V • u,t dx dt + f v(t*) • [u]‘+ dx = [ ([C/ltT 1 + R n ) dx . 

Jl n Jo Jo Jo 

Summing over all time slabs, the first term on the right-hand side of this equation 
vanishes except for initial and final time slab contributions. Utilizing nonnega- 
tivity of the remainder terms R n then yields the following inequality for the time 
evolution terms 

j v-u <t dxdt + j v(i") ■ [ u ]‘ t i dx^j > J (U(t”)-U(t°_)) dx. 
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Assume satisfaction of the system E-flux condition, the spatial term analysis used 
in the proof of Theorem 3.1 reduces to the inequality 

f ( f — v x ■ ■ U dx + f v_ • h ds'] dt > 0 . 
n=0 KeT K K ’ ’ JdK J 


Combining temporal and spatial results yields 

o = Bdg(v,v) > [ (U{t") - U{t°_ )) dx . 

Jo. 

Hence, the desired upper bound in (3.15) is established when applied to finite- 
dimensional subspaces 

(3.16) f U(u(v h (x,t^)))dx < [ I7(u(v/ l (x, )) dx . 

Jo Jo 

To obtain the lower bound in (3.15), we exploit the well-known thermodynamic 
concept of a minimum total entropy state (see for example [Mer88]). Define the 
integral average state u* at time slab boundaries 

U*(g) = — [ u(v h (x,tl))dx , n = 0,...,N. 
meas(ft) Jq 

For the DG space-time discretization of the Cauchy initial value problem, u* is 
invariant when evaluated at time slab boundaries, i.e. 


(3.17) u*(t”) = u’ttr 1 ) = . . . = u 


owing to discrete conservation in both space and time. A Taylor series with in- 
tegral remainder expansion of the entropy function given two states u*(£!l) and 
u(vh(x, t” )) for a fixed n yields 


U( u) = U(u*) +v(u*) • (u-u*) + [ (1 — 0))(u — u*) • (7 )UU (^)(u — u*) dO . 

Jo 

When integrated over ft, the second right-hand side term vanishes identically by 
the definition of u* 

[ U (u) dx = [ U{u*)dx + [ f \l -9)){u -u*) -U uli {9){u - u*) dOdx . 
JO JO JO Jo 

From strict convexity of the entropy function, it follows that u* is a minimum total 
entropy state since f Q U dx is minimized when u — u*. Finally, since u *(£* ) is 
constant for n = 0, . . . , N, then 

f U(u*(t°_))dx= [ U(u*(t?))dx< f U(u(v h {x, t?)))dx. 

Jo Jo Jo 

This establishes the lower bound in (3.15). □ 
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3.2. DG Stability Analysis for Systems with Soienoidal Involution. Our 

attention shifts to the MHD system with soienoidal involution 

f u t + ft ,* t = 0 
(3.18) < B i}Xi =0 

{ u(z,t°) = U 0 (x) 

with convex entropy extension 


(3.19) 


U } t + Fi iXi < 0 . 


The goal is to derive sufficient conditions for MHD system discretizations so that 
the cell entropy inequality (3.5), the global semi-discrete bound (3.6), and the 
global space-time bound (3.7) are obtained. Motivated by the Godunov MHD 
symmetrization theory, we consider an implementation of the DG method using 
the Godunov augmented MHD system. 

DG FEM for MHD: Find e V h such that 

(3.20) 5DG-MHD(Vh,W h ) =0 , Vw k e V- fe 

with 


SDOMHD(v 1 w)=vf £ [ /-(u(v) • w, t + fi(v) • W' Xi )dxdt 
n=0\K & T JlnJl < 


(3.21) 


- Y2 j f ok (w ■ <j> T v ) V • B(v)dxdt 
Ker J * njK 

-h / / w(a:_) • h(v(x_), v(x + );n ) ds dt 

Ker Jl- JdK 

+ Y, j (w(C +I )-u(v(C + 1 ))-w(^)-u(v(e))) dx\ . 


Observe the added V • B term with adjustable coefficient ok is motivated by the 
theory given in Sect. 2.2. The value of ok will be determined from the discrete 
energy analysis. This term is identical to that proposed by Powell [Pow94] using a 
different motivating argument. Unfortunately, without placing further constraints 
on the discrete B field, the Powell term is only valid for classical (smooth) solu- 
tions since this term cannot be written in divergence form. Consequently incorrect 
Rankine-Hugonoit jump conditions are observed for computed weak (discontinu- 
ous) solutions [Csi02]. Note that this term vanishes identically and correct weak 
solutions are obtained when a locally divergence- free basis is employed. 

A DG analysis similar to that used in Theorem 3.1 yields the following con- 
ditions for a discrete cell entropy inequality for the MHD formulation. 

Theorem 3.3 (DG Semi-Discrete MHD Cell Entropy Inequality). Let 
v/j 6 V c h denote a numerical solution obtained using the discontinuous Galerkin 
method (3.21) assuming continuous in time approximation for the MHD Cauchy 
initial value problem (3.18) with convex entropy extension (3.19). Assume the 
following conditions are satisfied: 
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1. Either a k = 1 or the local solenoidal condition holds pointwise 

V-B(v h )| K «0. 

2. The MHD system E-flux condition 

M± • (h(v_,v + ;n)-f(v(0))-n + <£(v(0)) (B(v(0)) • n)^) <0, 
V# € [0, 1] where v(0) = «f 0 [v] 

The numerical solution then satisfies the local semi-discrete cell entropy in- 
equality 

(3.22) f U(vh)dx+ [ F(v-,h, v + ft;n)ds < 0 , foreachKeT 

JK JdK 

with 

(3.23) F(v_,v + ;n) = (v)t • h(v_,v + ;n) - (F • n - <t> B • n)± 
as well as the global semi-discrete entropy inequality 

(3.24) j t jjJ{w h )dx< 0. 


Proof. Evaluate the energy, I?dg( v /i, v/J, for a single stationary element 
K in the DG discretization of the MHD system assuming continuous in time 
approximation 


A 

dt 


L 


U dx = 


— / (— v )iBi • fi)dz 4- / v_ • h ds 

— [ (— ■ n + <j)~ (B_ • n) + v • h) ds 
JdK 

— I (1 — ctk) <f> V ■ B dx 
Jk 

— [ (F(v_,v + ;n) + L>(v_,v + ;n))ds 
JdK 

— f (1 - cr/c) 4> V • B dx . 

Jk 


The remaining element interior term vanishes identically by either imposing 
crx — 1 or the local solenoidal condition on the magnetic induction field, 
V ■ B\ k = 0. Suitable definitions for the conservative entropy flux and entropy 
dissipation are given by 


.F(v_, v+; n) = (v)t * h(v_, v+; n) + {-T • n + <j> B ■ n)t 
£>(v_,v + ; n) = — 3([v]± • h + [-F ■ n + <f> B • n]+) . 

This choice of numerical entropy flux satisfies conservation and consistency prop- 
erties 
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• F(v_, v + ;n) = -F(v+, v_; -n) (conservation) 

• F(v, v; n) = (v ■ f — T + (j> B) • n = F • n (consistency) . 

Rewriting the jump term appearing in the entropy dissipation term as a path inte- 
gration assuming a parameterized state space v(0) — v_ +9 [v] 1 

£>(v_, v+; n) - -l([v]l • h + \-T • n + <j> B ■ n]t) 

= -|[v]i • (h - jf 1 (^(v(0)) ■ n - (</> B • n) T v (v(0))) dfl) 

M± ■ (h - f(v(<9)) ■ n + <f>{v{e)) (B(v(0)) • n)£) d9 . 

A sufficient condition for nonnegativity of D(v_, v+; n) is that the integrand be 
nonpositive. This yields the MHD E-flux condition 

[v]t-(h(v_,v + ;n)-f(v(0)).n+d>(v(<9)) (B(v{6))-n) T v ) < 0 , V0 € [0, 1] . 

This establishes the semi-discrete cell entropy inequality for MHD. Summation 
of (3.22) over all elements in the mesh together with the conservative telescop- 
ing property of F(v„, v+;n) yields the global semi-discrete entropy inequality 

(3.24) . □ 

The conditions set forth in Theorem 3.3 are also sufficient to establish two-sided 
bounds on the total entropy. 

THEOREM 3.4 (DG Fully-discrete MHD Total Entropy Bounds). Let 

Vh G V h denote the space- time numerical solution obtained using the discontin- 
uous Galerkin method (3.21) for the MHD Cauchy initial value problem (3.18) 
with convex entropy extension (3.19). Assume the following conditions are satis- 
fied: 

1. Either gk = 1 or the local solenoidal condition holds pointwise 

V-B(v fc )|* = 0. 

2. The MHD system E-flux condition 

[v]t • (h(v_, v + ; n) - f(v(0)) • n + <j>{v{6)) (B (v(0)) • n)£) < 0 , 

V0 G [0, 1] where v(0) — v_ + 6 [v]l. 

The numerical solution then satisfies the total entropy bound 

(3.25) f U(u*(t°_ )) dx < ( U{n{wh{x,t 1 l)))dx < f J7(u(v^(x, £°)) dx 

Jn Jn Ja 

where denotes the minimum total entropy state of the initial projected data 

= —2— l u (v h (x,t°_))dx . 
meas(iZ) Jq 



Proof. Omitted, see Theorem 3.2. □ 
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3.2.1. A Compatible B Field Representation. Unfortunately, conventional 
system E-fluxes do not satisfy the MHD system E-flux condition. Furthermore, 
calculation of the actual symmetrization variables for the MHD system (2.2) as- 
sociated with the entropy function, t/(u) — — ps, reveals that B is not a vector 
component of v, viz. 


(3.26) 


v=£/ u T = (7-1) 


/ 0=£ + £Xi ' 

I 7-1 + 2 p 
£V 

le. 

v 

eS. 

P 


Observe, however, that the last vector component pB/p is a — B multiple of the 
preceding component —p/p. Hence, it is possible to parameterize v on a line, 
v ($) = v_ “F 6 [v]l, and constrain B • n independent of 0 so that [B • n]i = 0. 
The following lemma states that under this constraint, the MHD system E-flux 
condition reduces to a constrained variant of the system E-flux condition (3.8). 

Lemma 3. 1 (B Field Compatibility). Assume the MHD system E-flux con- 
dition as given in Theorems 3.3 and 3.4. In addition , assume that B(v) • n is 
constrained to be continuous at interelement interfaces , i.e. [B(v) • n]t =0. 
Then, under this assumption , the results of Theorems 3.3 and 3.4 are identically 
obtained with the MHD system Eflux condition 


[v]t • (h - f(v(0)) • n + «£(v(0))(B(v(0)) • n) T v ) < 0 , V0e[0, 1] 


replaced by the constrained system Eflux condition 


[v]l • (h - f(v(0)) • n)| B •n const 0 , vee [0,i] . 


Proof. The result follows immediately since 

(3.27) Ml ■ (B(v(fl)) ■ n) T v = = 0 

due to the 0 independence of B • n at element interfaces. □ 

This result indicates the underlying intrinsic compatibility requirement of 
continuity in the normal component of the magnetic induction field for DG dis- 
cretizations of MHD. Precise implementational details are given in a separate 
work [Bar04]. In that same work, several other DG discretization formulations 
and simplified flux functions are given which satisfy the sufficient conditions 
given in Theorems 3.3 and 3.4 

• Transformed variable formulations 

• Constrained formulations 

• Penalty formulations 
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4. Conclusions. The energy analysis presented herein reveals the subtle in- 
terplay of involutions in the nonlinear stability of the DG method. Sufficient con- 
ditions for energy stability of DG discretizations of Maxwell and MHD systems 
have been obtained. From the viewpoint of discrete energy stability, analysis in- 
dicates that “standard” DG discretization Maxwell’s equations are energy stable 
without modification. Surprisingly, sufficient conditions for MHD discretization 
stability place more demanding requirements as set forth in Theorems 3.3 and 3.4. 
More complete details and DG formulations for MHD can be found in [Bar04]. 


REFERENCES 


[Bar98] 


[Bar99j 

[Bar04] 

[BB80] 

[BK04] 

[BoiSSj 

[Bos98] 

[BR02] 

[CHS90] 

[CLS89] 

[CLS04] 

[CS97] 

[Csi02] 

[Daf86j 

[DKK+02] 


T.J. Barth. Numerical methods for gasdynamic systems on unstructured meshes. In 
Kroner, Ohlberger, and Rohde, editors. An Introduction to Recent Developments 
in Theory and Numerics for Conservation Laws , volume 5 of Lecture Notes in 
Computational Science and Engineering, pages 195-285. Springer-Verlag, Hei- 
delberg, 1998. 

T.J. Barth. Simplified discontinuous Galerkin methods for systems of conservation 
laws with convex extension. In Cockbum, Kamiadakis, and Shu, editors, Discon- 
tinuous Galerkin Methods , volume 1 1 of Lecture Notes in Computational Science 
and Engineering. Springer-Verlag, Heidelberg, 1999. 

T.J. Barth. On the discontinuous Galerkin approximation of compressible ideal mag- 
netohydrodynamics I: Energy stable discretizations. In preparation , 2004. 

J.U. Brackbill and D.C. Barnes. The effect of nonzero V • B on the numerical solution 
of the magnetohydrodynamic equations. J. Comp. Phys., 35:426—430, 1980. 

N. Besse and D. Kroner. Convergence of locally divergence-free discontinuous 
Galerkin methods for the induction equations of the MHD system. Technical 
Report Submitted to M2AN, Wolfgang Pauli Institute, Austria, 2004. 

G. Boillat. Involutions des systems conservatifs. C. R. Acad. Sci. Paris, Sire /, 
307:891-894, 1988. 

A. Bossavit. Computational Electromagnetism, Variational Formulations, Comple- 

mentarity, Edge Elements. Academic Press, San Diego, 1998. 

P.B. Bochev and A.C. Robinson. Matching algorithms and physics: Exact sequences 
of finite element spaces. In D. Estep and S. Tavener, editors. Collected Lec- 
tures on the Preservation of Stability Under Discretization, Philadephia, 2002. 
SIAM. 

B. Cockburn, S. Hou, and C.W. Shu. TVB Runge-Kutta local projection discontinuous 

Galerkin finite element method for conservation laws IV: The multidimensional 
case. Math. Comp., 54:545-581, 1990. 

B. Cockbum, S.Y. Lin, and C.W. Shu. TVB Runge-Kutta local projection discontin- 
uous Galerkin finite element method for conservation laws III: One dimensional 
systems. J. Comp. Phys., 84:90-113, 1989. 

B. Cockburn, E Li, and C.W. Shu. Locally divergence- free discontinuous Galerkin 
methods for Maxwell equations. J. Comp. Phys., 194:588-610, 2004. 

B. Cockbum and C.W. Shu. The Runge-Kutta discontinuous Galerkin method for con- 

servation laws V: Multidimensional systems. Technical Report 201737, Institite 
for Computer Applications in Science and Engineering (ICASE), NASA Langley 
R.C., 1997. 

A. Csik. Upwind Residual Distribution Schemes for General Hyperbolic Conservation 
Laws with Application to Ideal Magnetohydrodynamics. PhD thesis. University 
of Leuven, Belgium, 2002. 

C. Dafermos. Quasilinear hyperbolic systems with involutions. Arch. Rational Mech. 

Anal, 106:373-389, 1986. 

A. Dedner, F. Kemm, D. Kroner, C.-D. Munz, T. Schnitner, and M. Wesenberg. Hyper- 



THE ROLE OF INVOLUTIONS IN DISCONTINUOUS G ALERKIN DISCRETIZATION 1 9 


[EMRS92] 

[God61] 

[God72] 

[HLvL83j 

[HS99] 

[JJS95] 

[JP86] 

[JS94] 

[LR74] 

[Mer88] 

[Moc80] 

[Ned80] 

[Osh84] 

[Pow94] 

[RH73] 

IShu99] 

(TOO) 

[Yee66] 


bolic divergence cleaning for the MHD equations. J. Comp. Phys 175:645-673, 
2002. 

B. Einfeidt, C. Munz, P. Roe, and B. Sjogreen. On Godunov -type methods near low 

densities. J . Comp . Phys., 92:273-295, 1992. 

S. K. Godunov. An interesting class of quasilinear systems. Dokl. Akad. Nauk. SSSR, 
139:521-523, 1961. 

S. K. Godunov. The symmetric form of magnetohydrodynamics equation. Num. Meth. 
Mech. Com. Media , 1:26-34, 1972. 

A. Harten, P. D. Lax, and B. van Leer. On upstream differencing and Godunov-type 
schemes for hyperbolic conservation laws. SIAM Rev., 25:35-61, 1983. 

J.M. Hyman and M. Shashkov. Mimetic discretizations for Maxwell’s equations. J. 
Comp. Phys., 151:881-909, 1999. 

J. Jaffre, C. Johnson, and A. Szepessy. Convergence of the discontinuous Galerkin 
finite element method for hyperbolic conservation laws. Math. Models and Meth- 
ods in Appl Set, 5(3):367-386, 1995. 

C. Johnson and J. Pitkaranta. An analysis of the discontinuous Galerkin method for a 

scalar hyperbolic equation. Math. Comp., 46:1-26, 1986. 

G. Jiang and C.-W. Shu. On a cell entropy inequality for discontinuous galerkin meth- 
ods. Math. Comp., 62:531-538, 1994. 

P. LeSaint and P.A. Raviart. On a finite element method for solving the neuton trans- 
port equation. In C. de Boor, editor, Mathematical Aspects of Finite Elements in 
Partial Differential Equations, pages 89-145. Academic Press, 1974. 

M. L. Merriam. An Entropy-Based Approach to Nonlinear Stability. PhD thesis, 
Stanford University, 1988. 

M. S. Mock. Systems of conservation laws of mixed type. J. Di iff. Eqns., 37:70-88, 
1980. 

J. C. Nedelec. Mixed finite elements in H 3 . Numer. Math., 35:315-341, 1980. 

S. Osher. Riemann solvers, the entropy condition, and difference approximations. 
SIAM J. Numer. Anal. , 2l(2):2i7-235, 1984. 

K. G. Powell. An approximate Riemann solver for magnetohydrodynamics (that works 

in more than one dimension). Technical Report 94-24, Institute for Computer 
Applications in Science and Engineering (ICASE), NASA Langley R.C., 1994. 
W. H. Reed and T. R. Hill. Triangular mesh methods for the neutron transport equation. 
Technical Report LA-UR-73-479, Los Alamos National Laboratory, Los Alamos, 
New Mexico, 1973. 

C.-W. Shu. Discontinuous Galerkin methods for convection-dominated problems. In 
Barth and Deconinck, editors, High-Order Discretization Methods in Computa- 
tional Physics , volume 9 of Lecture Notes in Computational Science and Engi- 
neering. Springer- Verlag, Heidelberg, 1999. 

G. Toth. The V * B = 0 constraint in shock-capturing magnetohydrodynamics codes. 
J. Comp. Phys., 161:605-652, 2000. 

K.S. Yee. Numerical solution of initial boundary value problems involving Maxwell’s 
equations in isotropic media. IEEE Trans. Ant. Prop., AP- 14:302-307, 1966. 


